Method And System For Error Compensation

ABSTRACT

A method for generating a set of kernels for convolution error compensation of a projection image of a physical object recorded by an imaging system comprises calculating the set of kernels in such a way that for each pixel of the projection image an asymmetric scatter distribution for error compensation is calculated representing a X-ray scatter originating along a ray from an X-ray source to the pixel.

The invention relates to a method for generating a set of kernels, amethod and a system for error compensation, a computer readable mediumand a program element, in particular to a method for convolution-basederror compensation of X-ray scatter.

Computed tomography (CT) is a process of using digital processing togenerate a three-dimensional image of the internal of an object underinvestigation (object of interest, object under examination) from aseries of two-dimensional x-ray images taken around a single axis ofrotation. The reconstruction of CT images can be done by applyingappropriate algorithms.

A basic principle of CT imaging is that projection data of an objectunder examination are taken by detectors of a CT system. The projectiondata represent information of the object passed by radiation beams. Togenerate an image out of the projection data these projection data (lineintegrals) can be back-projected leading to a two-dimensional image,i.e. representing a disc. Out of a plurality of such two-dimensionalimages a so-called voxel representation, i.e. a representation ofthree-dimensional pixels, can be reconstructed. In case that thedetectors are already arranged in form of a plane, two-dimensionalprojection data are achieved and the result of the back-projection is athree-dimensional voxel. That is, in modern, more sophisticatedso-called “cone-beam” CT and reconstruction methods the projection dataof two-dimensional detectors, i.e. detectors having a plurality ofdetecting elements arranged in form of a matrix, are directlyback-projected into a three-dimensional distribution of voxels in onesingle reconstruction step.

Scattered radiation is a major source of artifacts in cone-beam X-raycomputed tomography. By causing artifacts such as noise, streaks andlow-frequent inhomogeneities, so-called cupping, in reconstructedimages, scatter impedes visibility of soft contrasts, i.e. portionshaving low contrasts. Especially in volume imaging using interventionalX-ray systems where anti-scatter grids are inefficient, reliable andaccurate retrospective methods for scatter compensation are needed. Oneapproach for a correction are the so-called convolution based methodswhich are frequently used to estimate the scatter background inradiographic images. For example, such convolution based methods aredescribed in “Computerized scatter correction in diagnostic radiology”,K. P. Maher and J. F. Malone, Contemporary Physics 38(2), 131-148, 1997.

Although these convolution based correction methods do increase thequality of the reconstructed images, the reconstructed images stillexhibit artifacts, in particular in volumetric images.

It may be desirable to provide an alternative method for generating aset of kernels, a method and a system for error compensation, a computerreadable medium and a program element which may exhibit greater accuracyin error compensation or may be less prone to artifacts in thereconstructed image.

This need may be met by a method for generating a set of kernels, amethod and a system for error compensation, a computer readable mediumand a program element according to the independent claims.

According to an exemplary embodiment of a method for generating a set ofkernels for convolution error compensation of a projection image of aphysical object recorded by an imaging system, the method comprisescalculating the set of kernels in such a way that for each pixel of theprojection image a generally asymmetric scatter distribution for errorcompensation is calculated representing the X-ray scatter originating ina volume defined by the beam between an X-ray source and the pixel.

According to an exemplary embodiment a method for error compensation ofan image of a physical object comprises receiving an original projectionimage of an imaged physical object, converting the original projectionimage into a water-equivalent image (an image of water-equivalentthickness), extracting a number of scalar parameters from saidwater-equivalent image and possibly its gradient, determining at leastone pre-calculated kernel according to an exemplary embodiment of themethod for generating a set of kernels by relating the extractedparameters to the parameters of the pre-calculated kernels, andcompensating an error of the original projection image by using thedetermined at least one pre-calculated kernel.

According to an exemplary embodiment a system for error compensation ofan image of a physical object comprises a receiving unit adapted toreceive an original projection image of an imaged physical object, acalculation unit adapted to convert the original projection image into awater-equivalent image, in particular calculate the correspondinggradient image, and to extract a number of parameters from the images ofwater-equivalent thickness and in particular of the gradient, andetermination unit adapted to determine at least one pre-calculatedkernel according to an exemplary embodiment of the method for generatinga set of kernels by relating the extracted parameters to the parametersof the pre-calculated kernels, and a compensation unit adapted tocompensate an error of the original projection image by using thedetermined at least one pre-calculated kernel.

According to an exemplary embodiment a tomography apparatus comprises aradiation source, a radiation detector, and a system for errorcompensating according to an exemplary embodiment, wherein the radiationdetector is adapted to record data representing information of theoriginal projection image of the imaged physical object.

According to an exemplary embodiment a computer readable medium isprovided in which a program for error compensation of an image of aphysical object is stored, which program, when executed by a processor,is adapted to control a method comprising: receiving an originalprojection image of an imaged physical object, converting the originalprojection image into a water-equivalent image, determining at least onepre-calculated kernel according to an exemplary embodiment of the methodfor generating a set of kernels by relating the extracted parameters tothe parameters of the pre-calculated kernels, and compensating an errorof the original projection image by using the determined at least onepre-calculated kernel.

According to an exemplary embodiment a program element for errorcompensation of an image of a physical object is provided, whichprogram, when executed by a processor, is adapted to control a methodcomprising: receiving an original projection image of an imaged physicalobject, converting the original projection image into a water-equivalentimage in particular calculate the corresponding gradient image,extracting a number of scalar parameters from said water-equivalentimage and in particular from the gradient image, determining at leastone pre-calculated kernel according to an exemplary embodiment of themethod for generating a set of kernels by relating the extractedparameters to the parameters of the pre-calculated kernels, andcompensating an error of the original projection image by using thedetermined at least one pre-calculated kernel according to an exemplaryembodiment of the method for generating a set of kernels.

It may be seen as the gist of an exemplary embodiment of the presentinvention that a method for pre-calculating kernels is provided, whichkernels adequately accounts for the asymmetry of scatter distributionsgenerated along a ray, in dependence of the position where the raypenetrates the images object. One exemplary aspect of the presentinvention may be seen in that the present invention accurately accountsfor the fact that a large fraction of detected scattered X-ray quantamay originate from regions near the boundary of the imaged physicalobject and that the scatter distribution generated along the path tosuch locations may be highly asymmetric. The exemplary embodiment mayprovide a correction scheme that may offer the potential to much morequantitatively estimate and correct for scatter in radiographic imagesand projections of cone-beam computer tomography (CT) acquisitions.Thereby, possibly reducing image artifacts and thus possibly enhancinglow-contrast visibility, compared to a convolution-based method whichdoes not take into account a dependence on the position where the raypenetrates the images object, e.g. whether the considered pixel relatesto a centre of the physical object or to a border region. Preferably,the calculation of the set of kernels is done in such a way that foreach pixel of the projection image an asymmetric scatter distributionfor error compensation is calculated representing a X-ray scatteroriginating along a ray from an X-ray source to the pixel, whereinasymmetric may mean that no symmetry axis is existing. In particular,this asymmetry may be existing even in the case no anti-scatter grid isused.

The provided convolution-related scatter estimation scheme (not based onconvolution in the strict mathematical sense) uses pre-calculatedscatter kernels that determine the scatter contribution of a ray fromthe X-ray source to a detector element, depending on the objectattenuation at that pixel, and on further properties derived from theprojection image, such as estimates of the total object size or itsmaximal depth, or of the attenuation gradient in the water-equivalentimage at said pixel. The total scatter image may be obtained by summingup the contributions of all such rays. The kernels may be generatedeither experimentally or numerically. These kernels may be usable inorder to error compensating a projection image P⁽⁰⁾(x,y) which iscomprised of a primary portion P and a scatter portion S, i.e.P⁽⁰⁾(x,y)=P(x,y)+S(x,y).

This reconstruction method may be usable in the field of tomographyapparatuses, e.g. a computed tomography apparatus, in particular in anX-ray computer tomography.

In the following, further exemplary embodiments of the method forgenerating a set of kernels will be described. However, theseembodiments apply also for the method and the system for errorcompensation, the tomography apparatus, the computer readable medium andthe program element.

According to another exemplary embodiment of the method of generating aset of kernels the set of kernels is experimentally determined by usingan X-ray phantom as a model. In particular, in the calculating of theset of kernels results of an experimental measurement might be used.

According to another exemplary embodiment of the method for generating aset of kernels the set of kernels is calculated by using scattersimulations of a geometric model, preferably assuming water-likescattering characteristics, or scattering characteristics of othermaterials. Preferably, each kernel of the set of kernels is a functionof parameters of the geometric model.

That is, for the generation of pre-calculated kernels, normalizedscatter distributions K(x,y) may be off-line generated using pencil-beamMonte-Carlo scatter simulations of the geometric model, which may beparameterized in such a way that it takes into account a correct systemgeometry, e.g. geometry of a tomography system, a correct beam spectrum,e.g. the energy spectrum of the corresponding radiation source of thetomography system, and a correct anti-scatter grid, e.g. whether ananti-scatter grid and which specific anti-scatter grid is used in thetomography system. From these scatter distributions an estimation of ascatter image S⁰(x,y) may be obtainable by summing up the pre-calculatedcontributions for rays impinging on the individual detector pixels in aprojection image.

According to another exemplary embodiment of the method for generating aset of kernels at least one of the parameters is a radius of thegeometric model. Preferably, the kernel is further a function of a shiftbetween the projected centre of the geometric model and the positionwhere the penetrating pencil beam impinges onto the detector.

That is, for a given system configuration, e.g. tomography systemconfiguration, separate kernels K_(M,r,Φ)(x,y) may be off-line generatedas a function of model parameters M, e.g. at least one radius, and as afunction of a positional shift (r, Φ) of the model with respect to apencil beam used for the simulation, wherein (r, Φ) are polarcoordinates denoting the shift in a plane parallel to a detector planeof the tomography system. By calculating the kernels as a function of apositional shift (r, Φ) it may be possible to account for the scattervariation depending on the pixel location, e.g. whether the pixel is aboundary pixel or a centre pixel.

According to another exemplary embodiment of the method for generating aset of kernels the geometric model is an ellipsoidal model. Preferably,each kernel of the set of kernels is a function of r₁, r₂ and r₃ of thegeometric model and of a shift r, Φ between the centre of the model andthe position where the pencil beam penetrates the model possiblyresulting in a shift between the projected centre of the geometric modeland the position where the penetrating pencil beam impinges onto thedetector, wherein r₁, r₂ and r₃ may be the half axes of the ellipsoidalmodel.

For these model parameter M, e.g. r₁, r₂, and r₃ the pre-calculatedkernels may be calculated, i.e. the pre-calculated kernels K(x,y) may becalculated as a function of these model parameter M and as a function ofa positional shift (r, Φ) of the model with respect to the pencil beam.For each combination of model parameters M=(r₁, r₂, r₃), scatter kernelsK_(M,r,Φ)(x,y) may be generated under variation of the relative positionbetween pencil beam and ellipsoid model in the plane parallel to thedetector, wherein the positional shift of the model ellipsoid withrespect to the pencil beam is denoted by the polar coordinates (r, Φ).

According to another exemplary embodiment of the method for generating aset of kernels the geometric model is a spherical model. Preferably,each kernel of the set of kernels is a function of a radius R of thespherical model and of a shift r,Φ between the centre of the model andthe position where the pencil beam penetrates the model possiblyresulting in a shift between the projected centre of the geometric modeland the position where the penetrating pencil beam impinges onto thedetector.

In the following, further exemplary embodiments of the method for errorcompensation will be described. However, these embodiments apply alsofor the method for generating a set of kernels, the system for errorcompensation, the tomography apparatus, the computer readable medium andthe program element.

According to another exemplary embodiment of the method for errorcompensation the original projection image is normalized.

According to another exemplary embodiment of the method for errorcompensation each kernel of the set of kernels is a function of ageometry of the imaging system, a beam spectrum of the imaging systemand/or anti-scatter grid parameters of the imaging system.

In this context “normalized” denotes the fact that the quantity Pdenotes the detected intensity of primary radiation normalized by thevalue for air, so that P=1 corresponds to direct radiation and P=0corresponds to total absorption. By normalizing the projection image andconverting it to a water-equivalent image it may be possible to providefor an efficient way to error compensating images by using pre-calculatekernels.

According to another exemplary embodiment of the method of errorcompensation the original projection image is converted into awater-equivalent image according

${{T\left( {x,y} \right)} = \frac{- {\ln \left( {P^{(0)}\left( {x,y} \right)} \right)}}{\mu}},$

wherein P⁽⁰⁾ represents the original projection image, T(x,y) representsthe image of water-equivalent thickness T, and μ denoted the appropriateattenuation value of water.

In particular, generating a water-equivalent image may be suitable sincewater is predominant in a human being, thus leading to a simple butstill sufficiently correct model.

According to another exemplary embodiment the method for errorcompensation further comprises calculating a total scatter at a givenpixel of a pixel array by summing up the contributions of all kernelscorresponding to all pixels.

Such a summing up of all contributions of all kernels may be anefficient way to calculate the effect of the total scatteringcontribution on the measured intensity at the given pixel, e.g. theintensity measured at a detector element of a computer tomographyapparatus. This total scattering contribution may be afterwards used tocompensate errors introduced to the projection image by scattering.

According to another exemplary embodiment of the method for errorcompensation the total scatter at a given pixel is defined by:

${{S^{0}\left( {i,j} \right)} = {w \cdot {\sum\limits_{k,j}{K_{M,{r{({k,l})}},{\Phi {({k,l})}}}\left( {{i - k},{j - l}} \right)}}}},$

wherein: S⁰(i, j) is the total scatter at pixel (i,j), w denotes thearea of the pixel, and K_(M,r(k,l),Φ(k,l)) (i−k, j−l) is the kernelindicative for the scattering introduced from a ray impinging at pixel(k,l) at the location of pixel (i,j) and depending on: M, whichrepresents the parameters of a geometric model, and (r, Φ), whichrepresents a positional shift of the geometric model with respect to acentre of the pixel array.

According to another exemplary embodiment of the method for errorcompensation for the calculation of the kernels an ellipsoidal model isused, wherein M represents the half axes r₁, r₂, r₃ of the ellipsoidalmodel, wherein r₁=r₂=sqrt(A/π) and r₃=B, with A=a maximumcross-sectional area of the physical object, and with B=a maximumthickness of the physical object.

The corresponding parameters r₁, r₂ and r₃ may be extractable from theprojection image of the physical image and may represent the half axesof the ellipsoidal geometric model. Preferably, the projection image isconverted into an image of water-equivalent thickness T according

${{T\left( {x,y} \right)} = \frac{- {\ln \left( {P^{(0)}\left( {x,y} \right)} \right)}}{\mu}},$

wherein P⁽⁰⁾ represents the original projection image, wherein T(x,y)represents the image of water-equivalent thickness T; and μ denoted theappropriate attenuation value of water. The parameter A may be specifiedas the area of the shadow of the physical object on the projectionimage, e.g. the region in the projection with attenuation above acertain threshold, divided by the square of the system's geometricmagnification factor. The parameter B may be specified as the maximum ofT(x,y) after low-pass filtering or as a percentile from a histogram ofT, which both may minimize the influence of strong attenuationvariations. In an alternative embodiment, the model parameters may bedetermined from a least squares fit of a forward projection of the modelto the acquired projection.

It should be noted in this context that for a given model M, differentshift values are equivalent to different values of the water-equivalentthickness at the pencil beam position, ranging from the maximalthickness of the model at zero shift down to almost zero thickness atshifts almost equal to the spatial extent of the model. In turn, forsimple geometric models and a fixed shift angle Φ, a given value ofwater thickness T in the considered range unambiguously determines acorresponding value of r, so that in the interval (0,T_(max)], r(T) canbe assumed to be a unique function.

For this embodiment, consider the scatter contribution of a rayimpinging on the detector pixel with indices (k,l). At the location ofanother pixel (i,j), this ray produces a scatter contribution that isapproximately described by the expression K_(M,r(T(k,l)),Φ(k,l)) (i−k,j−l), where for the utilized kernel the shift radius r is specified bythe water thickness at pixel (k,l), and the shift angle Φ(k,l) might bechosen as the polar angle of pixel (k,l) in a coordinate system withorigin at the “centre of attenuation mass” (c₁,c₂) specified as

$\begin{pmatrix}c_{1} \\c_{2}\end{pmatrix} = {\frac{1}{\sum\limits_{k,l}{T\left( {k,l} \right)}} \cdot {\sum\limits_{k,l}{{T\left( {k,l} \right)} \cdot {\begin{pmatrix}k \\l\end{pmatrix}.}}}}$

The total scatter at pixel (i,j) may then be obtained by summing up thecontributions of all rays (k,l), yielding

${{S^{0}\left( {i,j} \right)} = {w \cdot {\sum\limits_{k,j}{K_{M,{r{({T{({k,l})}})}},{\Phi {({k,l})}}}\left( {{i - k},{j - l}} \right)}}}},$

where the sum runs over all pixel (k,l) of the detector, and w denotesthe pixel area.

According to another exemplary embodiment of the method for errorcompensation a spherical model is used for the calculation of thekernels, wherein the total scatter at the given pixel is defined by:

${{S^{0}\left( {i,j} \right)} = {w \cdot {\sum\limits_{k,j}{K_{{R{({{T{({k,l})}},{g{({k,l})}}})}},{r{({{T{({k,l})}}{g{({k,l})}}})}},{\Phi {({k,l})}}}\left( {{i - k},{j - l}} \right)}}}},$

wherein S⁰(i,j) is the total scatter at pixel (i,j), w denotes the areaof the pixel, and K_(R(T(k,l),g(k,l)),r(T(k,l)g(k,l)),Φ(k,l)) (i−k, j−l)is the kernel indicative for the scattering introduced from a rayimpinging at pixel (k,l) at the location of pixel (i,j) and dependingon: R, which represents a radius of the spherical geometric model, g,which represents a gradient of the corresponding image ofwater-equivalent thickness T, and (r, Φ), which represents a positionalshift of the ellipsoidal geometric model with respect to a centre of thepixel array.

According to another exemplary embodiment the method for errorcompensation method the parameters R and r are chosen according

$R = {\frac{T}{4} \cdot \sqrt{4 + g^{2}}}$

and

$r = {\frac{T}{4} \cdot g}$

and Φ=arg(grad T), with T=a water-equivalent thickness of the physicalobject, and g=|grad T|.

According to this exemplary embodiment a spherical geometric model maybe used, which may have a significant advantage in that it does notrequire estimation of global model parameters for each projection, butis based on local estimation of such parameters for each single ray.This variant uses spherical geometric models (phantoms) and, as thepreviously described variant, also works via phantom offsets withrespect to the pencil beam.

Applied to a projection P, the method first may require to calculate thegradient of the corresponding image of water-equivalent thickness T=−(lnP)/μ, which for each detector element exhibits a certain value ofmagnitude g=|grad T| and direction Φ=arg(grad T). To estimate thescatter contribution of a given source ray, the local values ofwater-equivalent thickness T, gradient magnitude g and direction Φ thenuniquely may determine the parameters (R, r, Φ) of the utilized spherephantom, where R may denote the radius of the sphere, and (r, Φ) may beits positional offset in the plane parallel to the detector. The mapping(T, g)

(R, r) is done in such a way that a parallel projection of the utilizedsphere would exhibit a water-equivalent thickness T and a thicknessgradient g at the position of the pencil beam. This is achieved by thefollowing equations:

${R = {\frac{T}{4} \cdot \sqrt{4 + g^{2}}}},\mspace{11mu} {r = {\frac{T}{4} \cdot {g.}}}$

It should be noted that in this way, the positional offset will be closeto zero in flat image areas, while it becomes close to the sphere radiusat steep gradients, e.g. near the object border. Using this method, fora given system geometry and beam quality, the convolution kernels arepre-calculated depending on the three parameters (R, r, Φ) as comparedto four parameters in case of the method based on ellipsoid models(kernels).

Such a spherical model may be in particular efficient when theprojection image is affected by truncations, e.g. in case the physicalobject is larger than the possible imaged object. While the ellipsoidalmodel may be affected by an erroneous estimation of the model parameterr₁=r₂ due to this truncations, the method based on sphere kernels may benot affected by truncations, due to its localized estimation of modelparameters.

According to another exemplary embodiment the method further comprisescalculating a first error compensated image in a multiplicative way byusing the total scatter. Preferably, the multiplicative correction isperformed according:

${P^{({n + 1})} = \frac{P^{(0)} \cdot P^{(n)}}{P^{(n)} + S^{(n)}}},$

wherein S^((n)) denotes the scatter image estimated from the projectionimage P^((n)).

The multiplicative way may in particular advantageous, since it mayexhibit an increased stability of convergence and may have theadditional advantage that negative projection values are avoided. Usingthe latter correction scheme, assuming the same estimated scatter, inregions with high attenuation a smaller amount of scatter may becorrected for as compared to regions with low attenuation. In contrastto subtractive correction where a pre-defined threshold on the maximalsubtracted amount of scatter may be specified in order to avoid negativeprojection values, such effects may be automatically avoided usingmultiplicative correction. In contrast to subtractive correction,multiplicative correction may need to be performed on full resolutionprojection images and therefore, in each iteration the estimated coarsescatter distributions may be again at least partly upsampled beforeapplying the correction step.

According to another exemplary embodiment the method further comprisescalculating a first error compensated image in a subtractive way byusing the total scatter. Preferably, the subtractive correction isperformed according: P^((n+1))=P⁽⁰⁾−S^((n)), wherein S^((n)) denotes thescatter image estimated from the projection image P^((n)).

According to another exemplary embodiment the method further comprisescalculating a second error compensated image by using the first errorcompensated image as the projection image. That is, the correction maybe performed in an iterative way, e.g. in 4 to 5 repetitions. That is,after estimation of a scatter image S⁽⁰⁾(x,y), this image is then usedto correct the originally acquired projection image P⁽⁰⁾(x,y)(containing both contributions of primary and scattered radiation),yielding an estimate P⁽¹⁾(x,y) of the true primary image. Because theinitial scatter-deteriorated projection image P⁽⁰⁾ results in a somewhatfalsified thickness image T, estimation and correction steps arepreferably repeated a number of times in an iterative fashion, untilconvergence of the estimated primary image is reached (this may usuallyachieved in about 4-5 iterations). Since scatter distributions aresmooth, scatter estimation may be carried out using a stronglydown-sampled detector pixel grid in order to decrease computationaleffort.

One exemplary aspect of the present invention may be seen in that avariable offset of the utilized phantoms is introduced during kernelgeneration. The schemes based on ellipsoid kernels and on sphere kernelsmake use of such an offset, and thus are potentially able toappropriately account for the asymmetry of scatter distributionsproduced near the object boundaries. Both estimation schemes may havehigh potential for application in X-ray volume imaging. Especially thescheme based on pre-calculated sphere kernels may produce accurateresults for different body regions (e.g., head, thorax and pelvisregions), and its performance may not be affected by the presence oftruncations. Most importantly, the optimal correction factors for thesebody regions may almost be the same. Regarding computational costs, thesphere method may be somewhat more demanding than the ellipsoidalmethods, since the scatter kernels of all possible sphere configurationsare preferably read and simultaneously stored in memory. For mostefficient use of this method, this data may be kept in memory instead ofrepeatedly being read when the method is applied to a sequence ofprojections of a rotational acquisition.

In order to improve the method based on ellipsoid kernels, which my beaffected by truncations when applied to the projections of thorax andpelvis, it may be possible to more robustly estimating the modelparameters via an optimization algorithm using forward projection of themodel, which might also be separately applicable to each acquiredprojection. This is due to the fact that this method relies on at leastapproximate estimation of two global parameters per projection, one ofwhich is difficult to estimate in case of truncations.

Furthermore, according to one exemplary aspect of the invention twodifferent schemes for the correction step of scatter compensation havebeen considered, namely subtractive and multiplicative correction. Eachscheme can be combined with each of the scatter estimation algorithmsaccording to exemplary embodiments. While subtractive correction maycasually produce clipping-related streak artifacts and may suffer fromunsatisfactory stability of the iterative estimation-correctionprocedure, it may be less computational time consuming. Alternatively,multiplicative correction may produce in all cases favourable results.Since multiplicative correction possibly needs to be performed on higherresolution projection images, in each iteration estimated coarse scatterdistributions may be again up-sampled before applying the correctionstep.

The error compensation of a projection image of a physical object may berealized by a computer program, i.e. by software, or by using one ormore special electronic optimization circuits, i.e. in hardware, or inhybrid form, i.e. by software components and hardware components. Thecomputer program may be written in any suitable programming language,such as, for example, C++ and may be stored on a computer-readablemedium, such as a CD-ROM. Also, the computer program may be availablefrom a network, such as the WorldWideWeb, from which it may bedownloaded into image processing units or processors, or any suitablecomputers.

It should be noted in this context, that the present invention is notlimited to computed tomography, but may include the use of C-arm based3D rotational X-ray imaging, positron emission tomography or the like.It should also be noted that this technique may in particular be usefulfor medical imaging of different body regions such as a head, a thoraxor a pelvic region of a patient.

These and other aspects of the present invention will become apparentfrom and elucidated with reference to the embodiment describedhereinafter. The disclosed embodiments and aspects described anywhere inthis application may be mixed and/or combined with each other.

An exemplary embodiment of the present invention will be described inthe following, with reference to the following drawings.

FIG. 1 shows a simplified schematic representation of a computedtomography system.

FIG. 2 shows a schematic sketch of a geometry for generation of anellipsoidal kernel.

FIG. 3 shows a schematic flow chart of an error compensation methodaccording to an exemplary embodiment of the invention.

FIG. 4 shows some exemplary scatter images.

The illustration in the drawings is schematically. In differentdrawings, similar or identical elements are provided with the similar oridentical reference signs.

FIG. 1 shows an exemplary embodiment of a computed tomography scannersystem which projection data may be handled by a correction methodaccording an embodiment of the invention.

The computed tomography apparatus 100 depicted in FIG. 1 is a cone-beamCT scanner. The CT scanner depicted in FIG. 1 comprises a gantry 101,which is rotatable around a rotational axis 102. The gantry 101 isdriven by means of a motor 103. Reference numeral 105 designates asource of radiation such as an X-ray source, which emits polychromaticor monochromatic radiation.

Reference numeral 106 designates an aperture system which forms theradiation beam emitted from the radiation source unit to a cone-shapedradiation beam 107. The cone-beam 107 is directed such that itpenetrates an object of interest 110 arranged in the center of thegantry 101, i.e. in an examination region of the CT scanner, andimpinges onto the detector 115 (detection unit). As may be taken fromFIG. 1, the detector 115 is arranged on the gantry 101 opposite to theradiation source unit 105, such that the surface of the detector 115 iscovered by the cone beam 107. The detector 115 depicted in FIG. 1comprises a plurality of detection elements 115 a each capable ofdetecting X-rays which have been scattered by, attenuated by or passedthrough the object of interest 110. The detector 115 schematically shownin FIG. 1 is a two-dimensional detector, i.e. the individual detectorelements are arranged in a plane, such detectors are used in so-calledcone-beam tomography.

During scanning the object of interest 110, the radiation source unit105, the aperture system 106 and the detector 115 are rotated along thegantry 101 in the direction indicated by an arrow 117. For rotation ofthe gantry 101 with the radiation source unit 105, the aperture system106 and the detector 115, the motor 103 is connected to a motor controlunit 120, which is connected to a control unit 125 (which might also bedenoted and used as a calculation, reconstruction or determinationunit).

In FIG. 1, the object of interest 110 is a human being which is disposedon an operation table 112. During the scan of a head 10 a, a thorax orany other part of the human being 110, while the gantry 101 rotatesaround the human being 110, the operation table 112 may displace thehuman being 110 along a direction parallel to the rotational axis 102 ofthe gantry 101. This may be done using a motor 113. By this, the head isscanned along a helical scan path. The operation table 112 may also bestopped during the scans to thereby measure signal slices.

The detector 115 is connected to the control unit 125. The control unit125 receives the detection result, i.e. the read-outs from the detectionelements 115 a of the detector 115 and determines a scanning result onthe basis of these read-outs. Furthermore, the control unit 125communicates with the motor control unit 120 in order to coordinate themovement of the gantry 101 with motors 103 and 113 with the operationtable 112.

The control unit 125 may be adapted for reconstructing an image fromread-outs of the detector 115. A reconstructed image generated by thecontrol unit 125 may be output to a display (not shown in FIG. 1) via aninterface.

The control unit 125 may be realized by a data processor or computer toprocess read-outs from the detector elements 115 a of the detector 115.

The computed tomography apparatus shown in FIG. 1 may capture computedtomography data of the head or thorax of a patient. In other words, whenthe gantry 101 rotates and when the operation table 112 is shiftedlinearly, then a helical scan is performed by the X-ray source 105 andthe detector 115 with respect to the patient. After having acquiredthese data, the data are transferred to the control unit 125, and themeasured data may be analyzed retrospectively.

FIG. 2 shows a schematic sketch of a geometry for generation of anellipsoidal kernel. By referring to this sketch an exemplary embodimentof ellipsoid kernels with variable offsets will be described. Thismethod accounts for the fact that the scatter contribution originatingfrom regions near the borders of the imaged object is highly asymmetric(e.g. object centre versus border region), while in known methods nooffset between scattering ray and model is used, and therefore theasymmetry of the produced scatter contribution is generally notaccurately accounted for.

According to the ellipsoidal model the projection image P(x,y) is firstnormalized and then converted into an equivalent image ofwater-equivalent thickness T(x,y) according to the formula T=−(ln P)/μ,where μ denotes the approximate attenuation value of water.

Afterwards, two scalars are extracted from the image of water-equivalentthickness T, specifying parameters of an ellipsoid-shaped model of theimaged object with water-like attenuation and scatteringcharacteristics. In particular, the homogeneous ellipsoid is assumed tohave half axes r₁=r₂=sqrt(A/π) in the plane parallel to the detectorsurface, and a half axis r₃=T_(max)/2 perpendicular to the detector.Here, A is a measure of the cross-sectional area of the imaged objectparallel to the detector surface and is specified as the area of theobject shadow (defined as the region in the projection withwater-equivalent thickness above a certain threshold, e.g. 10 mm)divided by the square of the system's geometric magnification factor.The quantity T_(max) is the approximate measure of largestwater-equivalent thickness. For calculation of scatter kernelsaccounting for the important dependence on the pixel position, knownwater slabs were replaced by the ellipsoid model, and additionallypositional offsets of the model with respect to the simulated pencilbeam were considered. For each combination of model parameters M=(r₁,r₂, r₃), scatter kernels K_(M,r,Φ)(x,y) were generated under variationof the relative position between pencil beam and ellipsoid model in theplane parallel to the detector, wherein the positional shift of themodel ellipsoid with respect to the pencil beam is denoted by the polarcoordinates (r, Φ). It should be noted in this context that for a givenmodel M, different shift values are equivalent to different values ofthe water-equivalent thickness at the pencil beam position, ranging fromthe maximal thickness of the model at zero shift down to almost zerothickness at shifts almost equal to the spatial extent of the model. Inturn, for a fixed shift angle f, a given value of water thickness T inthe considered range unambiguously determines a corresponding value ofr, so that in the interval (0,T_(max)], r(T) is a unique function. Thegeometry used for generation of ellipsoid kernels is illustrated in FIG.2.

FIG. 2 shows a flat detector 201 comprising columns x and rows y. On thedetector the detected scatter distribution of a ray penetrating anellipsoid model is schematically depicted by the white field on the flatdetector 201. It can be seen that the detected scatter is highlyasymmetric, i.e. the effect of scattering is much higher on the lefthalf of the flat detector 201 than on the right half of the flatdetector 201, leading to brighter pixels left from the centre. Further,FIG. 2 shows in a schematically shape the water ellipsoid 202 used togenerate the ellipsoid kernels K_(M,r,Φ)(x,y). The ellipsoid 202 ischaracterized by several parameters, in particular the half axes r₁ 203,r₃ 204, while r₂ is not depicted in FIG. 2 since it extendsperpendicular to the plane shown in FIG. 2. Furthermore a nonzeropositional shift r 205 is marked in FIG. 2, i.e. a nonzero shift of thefocal line 206, which extends from the focal spot 207 to the centre ofthe flat detector 201, and the centre of the water ellipsoid 202. As r₂the shift angle Φ is not shown in FIG. 2, since it is defined in a planeparallel to the detector surface. The water thickness T is depicted as208 in FIG. 2, while lines 209 schematically show different scatteredrays.

Using the model- and position-dependent kernels K_(M,r,Φ)(x,y), thescatter contribution of a ray impinging on the detector pixel withindices (k,l) at the location of another pixel (i,j) is given by theexpression K_(M,r(T(k,l)),Φ(k,l)) (i−k, j−l), where the length r of thepositional shift for the utilized kernel is specified via the waterthickness at pixel (k,l), and the shift angle Φ(k,l) is chosen as thepolar angle of pixel (k,l) in a coordinate system with origin at the“centre of attenuation mass” (c₁,c₂) specified as

$\begin{pmatrix}c_{1} \\c_{2}\end{pmatrix} = {\frac{1}{\sum\limits_{k,l}{T\left( {k,l} \right)}} \cdot {\sum\limits_{k,l}{{T\left( {k,l} \right)} \cdot {\begin{pmatrix}k \\l\end{pmatrix}.}}}}$

This may provide a suitable orientation of the asymmetric scatter kerneldistributions in case a single, approximately ellipsoid-shaped object.The total scatter at pixel (i,j) is then obtainable by summing up thecontributions of all rays (k,l) yielding the expression

${S^{0}\left( {i,j} \right)} = {w \cdot {\sum\limits_{k,j}{{K_{M,{r{({T{({k,l})}})}},{\Phi {({k,l})}}}\left( {{i - k},{j - l}} \right)}.}}}$

The geometry for sphere kernel generation is the same as the onedepicted in FIG. 2, except that r₁=r₂=r₃=R, i.e. instead of a waterellipsoid a water sphere is used. However, the offset may be calculateddifferently.

FIG. 3 shows a schematic flow chart of an error compensation methodaccording to an exemplary embodiment of the invention. This embodimentin particular relates to an ellipsoidal geometric model. The methodprocesses each acquired projection image separately and may comprise thefollowing sequence:

-   -   1. An acquired normalized projection image        P⁽⁰⁾(x,y)=P(x,y)+S(x,y) comprised of a primary portion P and a        scattered portion S is converted into an image of        water-equivalent thickness T according to T=−(ln P⁽⁰⁾ (x, y))/μ,        where μ denoted the appropriate attenuation value of water (Step        301).    -   2. From the image T, a number of scalar parameters are        extracted, specifying the parameters of a simple geometric model        of the imaged object. For instance, the half axes r_(i) of a        homogeneous ellipsoidal object model with circular cross-section        parallel to the detector plane and water-like attenuation may be        calculated according to, e.g. r₁=r₂=sqrt(A/π) and r₃=B, with        scalars A and B. In this particular embodiment, i.e. the        ellipsoid case, A is an approximate measure of the maximum        cross-sectional area of the imaged object parallel to the        detector surface, and B is an appropriate measure of the maximum        water-equivalent thickness of the imaged object. A may be        specified as the area of the object shadow, e.g. the region in        the projection with attenuation above a certain threshold,        divided by the square of the system's geometric magnification        factor. To minimize the influence of localized strong        attenuation variations, B may be specified as the maximum of        T(x,y) after low-pass filtering or as a percentile from a        histogram of T. In an alternative embodiment, the model        parameters are determined from a least square fit of a forward        projection of the model to the acquired projection (Step 302).    -   3. An estimation of the scattered image S⁽⁰⁾(x,y) is obtained by        summing up pre-calculated scatter contributions for the rays        impinging on the individual detector pixels. For this purpose,        normalized scatter distributions K(x,y) are off-line generated        using pencil-beam Monte-Carlo scatter simulations of the        parametric object model, taking into account the correct system        geometry, beam spectrum, and anti-scatter grid parameters. For a        given system configuration, separate kernels K_(M,r,Φ)(x,y) are        off-line generated as a function of model parameters M and, to        account for the important dependence on pixel position relative        to the projected object centre, as a function of a positional        shift (r, Φ) of the model with respect to the pencil beam,        wherein(r, Φ) are polar coordinates denoting the shift in the        plane parallel to the detector. It is important to note that for        a given model M, different shift values are equivalent to        different values of the water-equivalent thickness at the pencil        beam position, ranging from the maximal thickness of the model        at zero shift down to almost zero thickness at shifts almost        equal to the spatial extent of the model. In turn, for simple        geometric models and a fixed shift angle Φ, a given value of        water thickness T in the considered range unambiguously        determines a corresponding value of r, so that in the interval        (0,T_(max)], r(T) can be assumed to be a unique function.

Now, the scatter contribution of a ray impinging on the detector pixelwith indices (k,l) can be considered. At the location of another pixel(i,j), this ray produces a scatter contribution that is approximatelydescribed by the expression K_(M,r(T(k,l)),Φ(k,l)) (i−k, j−l), where forthe utilized kernel the shift radius r is specified by the waterthickness at pixel (k,l), and the shift angle Φ(k,l) might be chosen asthe polar angle of pixel (k,l) in a coordinate system with origin at the“centre of attenuation mass” (c₁,c₂) specified as

$\begin{pmatrix}c_{1} \\c_{2}\end{pmatrix} = {\frac{1}{\sum\limits_{k,l}{T\left( {k,l} \right)}} \cdot {\sum\limits_{k,l}{{T\left( {k,l} \right)} \cdot {\begin{pmatrix}k \\l\end{pmatrix}.}}}}$

The total scatter at pixel (i,j) is then obtained by summing up thecontributions of all rays (k,l), yielding

${{S^{0}\left( {i,j} \right)} = {w \cdot {\sum\limits_{k,j}{K_{M,{r{({T{({k,l})}})}},{\Phi {({k,l})}}}\left( {{i - k},{j - l}} \right)}}}},$

where the sum runs over all pixel (k,l) of the detector, and w denotesthe pixel area (Step 303).

-   -   4. Using the estimated scatter S⁰(x,y), the acquired projection        image P⁽⁰⁾(x,y) is then corrected, yielding an estimate        P⁽¹⁾(x,y) of the true primary image (Step 304). Because the        initial scatter-deteriorated projection image P⁽⁰⁾ results in a        somewhat falsified thickness image T, best results are achieved        when 1. to 4. (Steps 301 to 304) are repeated about four times        in an iterative fashion until convergence of the estimated        primary image is reached, wherein the repetition of 2. (Step        302) is optional. Since scatter distributions are smooth, the        above steps may be carried out using a strongly down-sampled        detector pixel grid in order to decrease computational effort.

Correction may either be performed in a subtractive or in amultiplicative way. Subtractive corrections in iteration n (n>1) iscarried out according to the formula P^((n+1))=P⁽⁰⁾−S^((n)). However,multiplicative correction according

$P^{({n + 1})} = \frac{P^{(0)} \cdot P^{(n)}}{P^{(n)} + S^{(n)}}$

was found to increase stability of convergence and has the additionaladvantage that negative projection values are avoided.

FIG. 4 shows some exemplary scatter images. In the upper figures of FIG.2 results of a error correction method according to an exemplaryembodiment of the present invention are shown in particular an ellipsoidmodel, whereas the lower figures of FIG. 2 show results of a knownmethod based on pre-calculated scatter kernels, which does not use apositional offset of the model and thus does not accurately account forthe asymmetric scatter contributions especially of the rays near theobject borders. In detail, FIG. 4 a shows in the upper part an estimatedscatter image depicted for a two-dimensional detector having rows y andcolumns x. In the lower part of FIG. 4 a the corresponding profile alongthe central horizontal cross-section through the image is displayed.FIG. 4 b shows in the upper part the estimated scatter image depictedfor a two-dimensional detector having rows y and columns x. In the lowerpart of FIG. 4 b the corresponding profile along the central horizontalcross-section through the image is displayed.

FIG. 4 c shows in the upper part a simulated ground truth depicted for atwo-dimensional detector having rows y and columns x. In the lower partof FIG. 4 a the corresponding profile along the central horizontalcross-section through the image is displayed. FIG. 4 d depicts the samefor the known method. FIGS. 4 c and 4 d are the same since the sameground truth is used for the comparison.

FIG. 4 e shows in the upper part the ration of estimated image to groundtruth in a two-dimensional plot having rows y and columns x. It can beclearly seen that the mean ratio is about 1 but still showing a slightlyoverestimation of about 5% while the scatter shape is well approximated,which can be seen from the relatively uniform distribution of thegreyscale values. On contrary, FIG. 4 f shows a highly non-uniformlydistribution. In particularin the image centre, corresponding to thearea of greatest thickness of the object, the scattering is highlyoverestimated, while the overestimation is much lower near the borders.In mean the overestimation is about 44%.

In the test implementation, the computational effort of the correctionmethod according to an exemplary embodiment of the invention was onlyabout twice as high than for the known convolution-based approaches, theresults of which are depicted in the lower part of FIG. 4. In generalthe correction method according to an exemplary embodiment of theinvention may allow for potentially much more accurate estimation of thescatter distribution, due to the fact that the dependence of scatteringon further parameters than only the water-equivalent thickness isaccounted for. It should be noted that the term “comprising” does notexclude other elements or steps and the “a” or “an” does not exclude aplurality. Also elements described in association with differentembodiments may be combined. It should also be noted that referencesigns in the claims shall not be construed as limiting the scope of theclaims.

1. A method for generating a set of kernels for convolution errorcompensation of a projection image of a physical object recorded by animaging system, the method comprising: calculating the set of kernels insuch a way that for each pixel of the projection image an asymmetricscatter distribution for error compensation is calculated representing aX-ray scatter originating in a volume defined by the beam between anX-ray source and the pixel.
 2. The method according claim 1, wherein theset of kernels is experimentally determined by using an X-ray phantom asa model.
 3. The method according to claim 1, wherein the set of kernelsis calculated by using scatter simulations of a geometric model.
 4. Themethod according to claim 2, wherein each kernel of the set of kernelsis a function of parameters of the geometric model.
 5. The methodaccording to claim 4, wherein at least one of the parameters is a radiusof the geometric model.
 6. The method according to claim 4, wherein thekernel is further a function of a shift between the projected centre ofthe geometric model and the position where the penetrating pencil beamimpinges onto the detector.
 7. The method according to claim 2, whereinthe geometric model is an ellipsoidal model.
 8. The method according toclaim 7, wherein each kernel of the set of kernels is a function of r1,r2 and r3 of the geometric model, and of a shift r,Φ between the centreof the model and the position where the pencil beam penetrates themodel.
 9. The method according to claim 2, wherein the geometric modelis a spherical model.
 10. The method according to claim 9, wherein eachkernel of the set of kernels is a function of a radius R of thespherical model and a shift r,Φ between the centre of the model and theposition where the pencil beam penetrates the model.
 11. The methodaccording to claim 1, wherein each kernel of the set of kernels is afunction of a geometry of the imaging system, a beam spectrum of theimaging system and/or anti-scatter grid parameters of the imagingsystem.
 12. A method for error compensation of an image of a physicalobject, the method comprising: receiving an original projection image ofan imaged physical object; converting the original projection image intoa water-equivalent image, in particular calculating the correspondinggradient image; extracting a number of parameters from the images ofwater-equivalent thickness and in particular from the gradient image;determining at least one pre-calculated kernel according to claim 1 byrelating the extracted parameters to the parameters of thepre-calculated kernels; and compensating an error of the originalprojection image by using the determined at least one pre-determinedkernel.
 13. The method according claim 12, wherein the originalprojection image is normalized.
 14. The method according claim 12,wherein the original projection image is converted into awater-equivalent image according${{T\left( {x,y} \right)} = \frac{- {\ln \left( {P^{(0)}\left( {x,y} \right)} \right)}}{\mu}},$wherein P(0) represents the original projection image, T(x,y) representsthe image of water-equivalent thickness T; and μ denotes the appropriateattenuation value of water.
 15. The method according to claim 12,further comprising: calculating a total scatter at a given pixel of anpixel array by summing up the contribution of all kernels correspondingto all pixels.
 16. The method according to claim 15 wherein the totalscatter at the given pixel is defined by:${{S^{0}\left( {i,j} \right)} = {w \cdot {\sum\limits_{k,j}{K_{M,{r{({T{({k,l})}})}},{\Phi {({k,l})}}}\left( {{i - k},{j - l}} \right)}}}},$ wherein: S⁰(i,j) is the total scatter at pixel (i,j), w denotes thearea of a pixel, and K_(M,r(T(k,l)),Φ(k,l)) (i−k, j−l) is the kernelindicative for the scattering introduced from a ray impinging at pixel(k,l) at the location of pixel (i,j) and depending on: M whichrepresents the parameters of the geometric model; and (r,Φ) whichrepresents a positional shift of the ellipsoidal geometric model withrespect to a centre of the pixel array.
 17. The method according toclaim 16 wherein for the calculation of the kernels an ellipsoidal modelis used; and wherein M represents the half axes r1, r2, r3 of theellipsoidal model wherein r1=r2=sqrt(A/π) and r3=B, with A=a maximumcross-sectional area of the physical object, and B=a maximum thicknessof the physical object.
 18. The method according to claim 15 wherein forthe calculation of the kernels a spherical model is used; and whereinthe total scatter at the given pixel is defined by:${{S^{0}\left( {i,j} \right)} = {w \cdot {\sum\limits_{k,j}{K_{{R{({{T{({k,l})}},{g{({k,l})}}})}},{r{({{T{({k,l})}}{g{({k,l})}}})}},{\Phi {({k,l})}}}\left( {{i - k},{j - l}} \right)}}}},$ wherein: S⁰(i,j) is the total scatter at pixel (i,j), w denotes thearea of the pixel, and K_(R(T(k,l),g(k,l)),r(T(k,l)g(k,l)),Φ(k,l)) (i−k,j−l) is the kernel indicative for the scattering introduced from a rayimpinging at pixel (k,l) at the location of pixel (i,j) and dependingon: R which represents a radius of the spherical geometric model; gwhich represents a gradient of the corresponding image ofwater-equivalent thickness T, and (r,Φ) which represents a positionalshift of the ellipsoidal geometric model with respect to a centre of thepixel array.
 19. The method of claim 18, wherein$R = {{{\frac{T}{4} \cdot \sqrt{4 + g^{2}}}\mspace{14mu} {and}\mspace{14mu} r} = {\frac{T}{4} \cdot g}}$ and Φ=arg(grad T), with T=a water-equivalent thickness of the physicalobject, and g=|grad T|.
 20. The method according to claim 16, furthercomprising: calculating a first error compensated image in amultiplicative way by using the total scatter.
 21. The method accordingclaim 20, further comprising: performing the multiplicative correctionaccording${P^{({n + 1})} = \frac{P^{(0)} \cdot P^{(n)}}{P^{(n)} + S^{(n)}}},$ wherein S(n) denotes the scatter image estimated from the projectionimage P(n).
 22. The method according to claim 16, further comprising:calculating a first error compensated image in a subtractive way byusing the total scatter.
 23. The method according claim 22, furthercomprising: performing the multiplicative correction according:P^((n+1))=P⁽⁰⁾−S^((n)), wherein S(n) denotes the scatter image estimatedfrom the projection image P(n).
 24. The method according to claim 20,further comprising: calculating a second error compensated image byusing the first error compensated image as the projection image.
 25. Asystem for error compensation of an image of a physical object, thesystem comprising: a receiving unit adapted to receive an originalprojection image of an imaged physical object; a calculation unitadapted to convert the original projection image into a water-equivalentimage, in particular to calculate the corresponding gradient image, andto extract a number of parameters from the images of water-equivalentthickness and in particular from the gradient image; a determinationunit adapted to determine at least one pre-calculated kernel accordingto claim 1 by relating the extracted parameters to the parameters of thepre-calculated kernels; and a compensation unit adapted to compensate anerror of the original projection image by using the determined at leastone pre-calculated kernel.
 26. A tomography apparatus comprising: aradiation source; a radiation detector; and a system for errorcompensating according claim 25; wherein the radiation detector isadapted to record data representing the original projection image of theimaged physical object.
 27. A computer readable medium in which aprogram for error compensation of an image of a physical object isstored, which program, when executed by a processor, is adapted tocontrol a method comprising: receiving an original projection image ofan imaged physical object; converting the original projection image intoa water-equivalent image, in particular calculating the correspondinggradient image; extracting a number of parameters from the images ofwater-equivalent thickness and in particular from the gradient image;determining at least one pre-calculated kernel according to claim 1 byrelating the extracted parameters to the parameters of thepre-calculated kernels; and compensating an error of the originalprojection image by using the determined at least one pre-calculatedkernel.
 27. A program element for error compensation of an image of aphysical object, which program, when executed by a processor, is adaptedto control a method comprising: receiving an original projection imageof an imaged physical object; converting the original projection imageinto a water-equivalent image, in particular calculating thecorresponding gradient image; extracting a number of parameters from theimages of water-equivalent thickness and in particular from the gradientimage; determining at least one pre-calculated kernel according to claim1 by relating the extracted parameters to the parameters of thepre-calculated kernels; and compensating an error of the originalprojection image by using the determined at least one pre-calculatedkernel.